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Abstract.  This  paper  reports  on  an  interdisciplinary  effort,  which  involves 
applied  mathematicians,  material  scientists  and  physicists  at  North  Carolina 
State  University,  to  integrate  new  intelligent  processing  approaches  with  ad¬ 
vanced  mathematical  modeling,  optimization,  and  control  theory  to  guide  the 
construction  and  experimental  implementation  of  a  series  of  high  pressure  (up 
to  100  atm)  organometallic  chemical  vapor  deposition  (CVD)  reactors.  An  in¬ 
tegral  component  of  this  research  program  is  the  design  of  the  reactor  so  that 
control  and  sensing  are  a  basic  component  of  the  optimal  design  efforts  for  the 
reactor.  We  report  here  on  the  successful  use  of  mathematics  in  a  fundamental 
role  in  the  development  of  linear  and  nonlinear  feedback  control  methods  for 
real-time  implementation  on  the  reactor.  This  is  achieved  in  the  required  con¬ 
text  of  gas  dynamics  coupled  with  nonlinear  surface  deposition  processes.  The 
problems  are  optimal  tracking  problems  (for  the  chemical  component  fluxes 
over  the  substrate)  that  employ  state-dependent  Riccati  gains  with  nonlin¬ 
ear  observations  and  the  resulting  dual  state  dependent  Riccati  equations  for 
the  compensator  gains.  This  control  methodology  is  successfully  combined 
with  reduced  order  model  methods  based  on  proper  orthogonal  decomposi¬ 
tion  techniques.  Computational  results  to  support  the  efficacy  of  our  approach 
and  methods  are  also  included. 


1.  Introduction 

Chemical  vapor  deposition  (CVD)  is  an  important  industrial  technique  used  to 
grow  thin  films  with  certain  desired  properties.  This  process  involves  the  deposi¬ 
tion  of  precursor  vapor  sources  onto  a  heated  substrate  where  they  react  to  form 
the  desired  material.  CVD  is  a  key  element  in  a  wide  variety  of  advanced  industrial 
applications,  ranging  from  the  development  of  short  wavelength  light  sources  to  de¬ 
tectors  and  integrated  sensors,  in  particular  the  integration  of  III-V  optoelectronics 
and  silicon  technology.  In  addition,  wide  bandgap  materials  are  of  particular  in¬ 
terest  for  advanced  silicon  ULSI  technology  in  the  context  of  dielectric  isolation, 
vertical  integration,  optically  interconnected  common  memory,  integrated  sensors 
and  microwave  applications.  These  materials  are  also  of  considerable  interest  in 
the  context  of  optoelectronics  for  applications  in  displays,  optical  recording,  signal 
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processing,  printing  and  medicine  [1].  In  these  processes,  the  control  of  the  con¬ 
centrations  and  distributions  of  both  point  defects  and  extended  defects  and  the 
related  control  of  the  surface  morphology  and  interfacial  chemistry  for  compound 
semiconductor  heterostructures  are  essential  because  of  their  important  role  in 
the  control  of  the  electrical/optical  properties  and  the  reliability  of  wide  bandgap 
semiconductor  devices  and  circuits. 

For  the  past  five  years,  an  interdisciplinary  effort  at  North  Carolina  State 
University  has  been  carried  out  to  explore  new  intelligent  processing  approaches 
that  access  conditions  outside  the  capabilities  of  conventional  methods.  Gener¬ 
ally  the  control  of  the  stoichiometry,  and  the  related  issue  of  the  control  of  the 
point  defect  chemistry  of  the  heteroepitaxial  layers  of  mixed  III-V  compounds 
requires,  under  the  conditions  of  thermally  activated  growth,  high  ratios  in  the 
fluxes  of  the  group  V  to  the  group  III  source  materials.  In  this  context  we  have 
explored  organometallic  CVD  (OMCVD)  under  superatmospheric  conditions  (up 
to  100  atm),  where  high  partial  pressure  ratios  can  be  established  without  com¬ 
promise  with  respect  to  the  growth  rate.  Because  the  process  is  operating  at  high 
flows/ vapor  densities,  control  of  the  fluid  dynamics  becomes  essential  for  optimal 
growth  conditions.  Therefore  advanced  methods  of  mathematical  modeling,  op¬ 
timization,  and  control  theory  have  been  applied  to  guide  the  development  and 
experimental  implementation  of  these  processes.  In  particular,  advanced  simula¬ 
tions  for  flow  processes  in  a  computer-aided  design  (CAD)  mode  resulted  in  several 
generations  of  reactor  geometries  before  a  suitable  configuration  that  promised  de¬ 
sirable  flow  characteristics  near  the  substrate  was  obtained.  In  this  presentation 
we  discuss  a  third  generation  reactor  design  (see  Fig.  1)  resulting  from  this  devel¬ 
opmental  process.  For  a  fourth  generation  reactor  design  see  [23] . 

The  mathematical  models  for  OMCVD  (under  superatmospheric  pressure) 
possess  one  of  the  most  complex  fluid  dynamics  systems  imaginable.  Some  of  the 
complex  issues  in  computing  chemically  reacting  flows  include  the  simulation  of 
three-dimensional  flows  governed  by  Navier-Stokes  equation  coupled  with  equa¬ 
tions  for  energy  and  species  and  strong  temperature  dependence  of  the  physical 
properties  of  gases.  Chemical  reactions  are  taking  place  in  the  gas-phase  as  well 
as  on  the  substrate.  However,  since  only  a  trace  amount  of  reactants  mixed  with 
carrier  gas  is  used,  a  dilute  approximation  is  assumed.  This  leads  to  a  quasi-steady 
gas-phase  model  with  steady-state  nonlinear  coupled  system  of  equations  for  the 
continuity,  momentum  and  energy  that  is  decoupled  from  the  time  dependent 
species  equations.  This  gas  phase  model  is  coupled  with  a  reduced  order  model  of 
the  surface  reactions  involved  in  the  decomposition  of  source  vapors  from  the  gas 
phase  and  the  growing  film  on  the  substrate. 

The  resulting  mathematical  model  is  a  system  of  partial  differential  equa¬ 
tions  coupled  with  nonlinear  ordinary  differential  equations  describing  the  surface 
deposition  process.  Numerical  simulations  and  control  designs  and  syntheses  of 
such  systems  are  faced  with  considerable  challenges  regarding  dimensionality  and 
nonlinearity.  This  paper  describes  our  efforts  during  the  last  five  years  to  overcome 
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these  difficulties.  More  specifically,  §2  describes  the  proper  orthogonal  decomposi¬ 
tion  (POD)  technique,  also  known  as  the  Karhunen-Loeve  procedure,  that  is  used 
to  obtain  low  dimensional  dynamic  models  of  distributed  parameter  systems.  The 
POD  method,  which  is  well  known  in  statistical  and  pattern  recognition  fields 
[2],  has  been  shown  to  be  an  effective  tool  for  the  analysis  of  complex  systems 
such  as  turbulence  flows,  shear  flows,  and  weather  prediction  (see  e.g.,  [3]  and  the 
references  therein).  Roughly  speaking,  POD  is  an  optimal  technique  of  finding  a 
basis  that  spans  an  ensemble  of  data,  collected  from  an  experiment  or  a  numerical 
simulation  of  a  dynamical  system,  in  the  sense  that  when  these  basis  functions 
are  used  in  a  Galerkin  procedure,  they  will  yield  a  finite  dimensional  system  with 
the  smallest  possible  degrees  of  freedom.  Thus  this  method  may  well  be  suited 
to  treat  optimal  control  and  parameter  estimation  of  distributed  parameter  sys¬ 
tems.  In  §3,  we  describe  our  successful  use  of  POD  techniques  as  a  reduced  basis 
method  for  computation  of  feedback  controls  and  compensators  in  a  high  pres¬ 
sure  CVD  reactor.  More  specifically,  we  present  a  proof-of-concept  computational 
implementation  of  this  method  with  a  simplified  growth  example  of  group  III-V 
compounds  that  includes  multiple  species  and  controls,  gas  phase  reactions  (no 
surface  reactions),  and  time  dependent  tracking  signals  that  are  consistent  with 
pulsed  vapor  reactant  inputs.  In  §4  state  estimation  and  feedback  tracking  con¬ 
trol  methods  for  nonlinear  systems  are  presented.  The  methods,  which  are  based 
on  the  “state-dependent  Riccati  equations”,  allow  the  construction  of  nonlinear 
estimators  and  nonlinear  feedback  tracking  controls  for  a  wide  class  of  systems 
including  high  pressure  CVD  systems  considered  here.  The  performance  of  the 
nonlinear  estimator  and  tracking  control  will  be  presented  on  a  flight  dynamics 
simulation  example.  Finally,  §5  contains  our  overall  conclusions. 


2.  Proper  Orthogonal  Decomposition 

In  general,  the  discretization  of  linear/nonlinear  partial  differential  equations  using 
finite  element,  finite  volume,  or  finite  different  methods  involves  basis  functions 
that  have  little  to  do  with  the  differential  equation.  For  example,  piecewise  poly¬ 
nomials  are  used  in  the  finite  element  method,  grid  functions  are  used  in  the  finite 
difference  method,  and  Legendre  or  Chebyshev  polynomials  are  used  in  some  spec¬ 
tral  methods.  POD,  on  the  other  hand,  uses  basis  functions  that  span  a  data  set, 
collected  from  an  experiment  or  numerical  simulation  of  a  dynamical  system,  in  a 
certain  “optimal”  fashion.  Because  POD  basis  elements  are  optimal  in  the  sense 
that  they  are  the  extractions  of  characteristic  features  of  the  data  set,  frequently 
only  a  small  number  of  POD  basis  functions  are  needed  to  describe  the  solution. 
POD  based  approximation  methods  have  been  applied  to  numerous  applications 
including  turbulent  coherent  flows  [4] ,  shear  flows  [5] ,  characterization  of  human 
faces  [6],  and  image  recognition  [7].  More  recently,  the  possibility  of  POD  based 
control  design  and  parameter  estimation  has  been  proposed.  In  particular,  appli¬ 
cations  of  POD  to  optimization  or  open  loop  control  were  developed  in  [8,  9,  10], 
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to  feedback  control  design  were  reported  in  [11,  12,  13,  14,  15],  and  to  parameter 
estimation  or  inverse  problems  were  discussed  in  [16].  References  to  recent  work 
of  other  authors  can  be  found  in  [11]- [16]  and  [23]. 

We  now  outline  an  algorithm  to  obtain  the  POD  basis  element.  The  mathe¬ 
matical  basis  for  the  algorithm  has  been  described  in  numerous  articles  (see  e.g., 
[8]).  Let  TO  f)  :  1  <  j  <  £  n]  denote  the  set  of  N  observations  (also 

called  snapshots)  of  some  physical  processes  over  a  domain  fl.  In  the  context  of 
CVD  process,  these  observations  could  be  experimental  measurements  or  numeri¬ 
cal  solutions  of  velocity  fields,  temperatures,  species  etc.  taken  at  different  physical 
parameters  (Reynolds  number,  input  flow  rates  etc.)  or  time  steps. 

Step  1.  Compute  the  covariant  matrix  C.  The  matrix  elements  of  C  are  given  by 

C'ifc  —  j  LJi(x)TJ^(^)d^r, 
for  i,  k  =  1,  2, . . . ,  N. 

Step  2.  Solve  the  eigenvalue  problem  CV  =  AV.  Since  C  is  a  nonnegative,  Her- 
mitian  matrix,  it  has  a  complete  set  of  orthogonal  eigenvectors 
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with  the  corresponding  eigenvalues  arranged  in  ascending  order  as  Ai  > 

A2  >  •  •  •  >  A n  —  0. 

Step  3.  Compute  the  POD  basis  vectors.  The  POD  basis  elements  4>,(a;)  such  that 
ArpOD  =  span{$i,  <I>2, . . . ,  dw},  where  XPOD  is  the  finite-dimensional 
POD  space,  are  given  as 

N 

$fc  =  5>iUi, 
i= 1 

where  1  <  k  <  N  and  a\  are  the  elements  of  the  eigenvector  Vfc  corre¬ 
sponding  to  the  eigenvalue  A*,. 

To  approximate  a  distributed  parameter  system  by  a  finite-dimensional  prob¬ 
lem  one  uses  a  combination  of  Galerkin  procedures  and  POD  basis  elements  (see 
[11,  12]  for  details).  However,  to  this  point  we  have  not  discussed  any  model  reduc¬ 
tion  features  associated  with  using  POD  basis  elements  in  approximation  schemes. 
In  the  algorithm  described  above,  the  number  N  may  be  large,  100  —  1000  or  even 
more,  depending  on  the  complexity  of  the  dynamics  represented  in  the  “snap¬ 
shots”  IR.  In  general,  one  should  take  N  sufficiently  large  so  that  the  snapshots 
Uj  contain  all  salient  features  of  the  dynamics  being  investigated.  Thus,  the  POD 
basis  functions  <!>,,  used  with  the  original  dynamics  in  a  Galerkin  procedure,  of¬ 
fers  the  possibilities  of  achieving  a  high  fidelity  model,  albeit  with  perhaps  a  large 
dimension  N. 
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To  achieve  model  reduction,  one  chooses  M  <C  N  and  carries  out  a  Galerkin 
procedure  with  the  set  of  elements  {$1,  <f>2, . . . ,  i\w}-  The  crucial  question  is  how 
to  choose  M.  As  discussed  elsewhere  (see  e.g.,  [11,  12])  the  percentage  of  the  total 
snapshots  set  data  variability  contained  in  a  certain  POD  mode  <&*,  is  given  by  the 
ratio  of  the  eigenvalue  A k  to  the  total  of  all  eigenvalues,  Afc/^yli  A j.  The  reason 
for  ordering  the  POD  modes  from  highest  to  lowest  eigenvalues  is  to  include  as 
much  of  the  variability  of  the  system  into  the  first  few  modes  as  possible.  Therefore 
to  capture  most  of  the  data  variability  of  the  system  contained  in  the  N  POD 
elements,  it  suffices  to  choose  M,  where  M  is  sufficiently  smaller  than  N,  so  that 
1  A  ~  Efci  A  j.  For  the  CVD  examples  studied  in  [8,  11,  12],  the  POD  system 
was  constructed  for  N  =  100  —  200  and  a  reduced  order  model  with  M  =  2  —  10 
resulted  in  a  truly  significant  computational  savings. 

While  the  above  comments  suggest  the  proper  choice  for  the  dimension  of 
the  reduced  order  model  to  be  used  in  simulations,  there  are  additional  order 
questions  related  to  the  linear  control  system  to  be  used  in  determining  reduced 
order  gains  and  compensators.  We  have  found  that  the  ranks  of  the  controllability 
and  observability  matrices  have  sometimes  been  useful  criteria  (see  the  discussion 
in  [13,  11,  12])  to  help  in  the  choice  of  the  number  of  modes  to  use  in  control 
design  applications  for  the  reduced  basis  representation. 

In  the  next  section,  we  present  a  proof-of-concept  computational  implemen¬ 
tation  of  this  method  with  a  simplified  growth  example  for  III-V  layers.  In  this 
example  we  implement  Dirichlet  boundary  control  of  dilute  reactants  transported 
by  convection  and  diffusion  to  an  absorbing  substrate  after  they  undergo  gas  phase 
reactions. 


3.  Application  of  POD  to  Compensator  Control  of  CVD 

The  particular  geometry  of  the  differentially  pressure  controlled  (DPC)  reactor 
system  under  consideration  here  features  horizontal  flow  of  the  process  gases  and 
source  vapor/carrier  gas  mixtures  into  an  expansion  section  leading  into  a  rect¬ 
angular  channel  that  contains  the  substrate  (see  Fig.  1).  The  substrate  wafer  is 
mounted  on  a  rotating  induction  heated  SiC  coated  graphite  susceptor.  The  ex¬ 
haust  gases  are  vented  through  a  vertical  exhaust  tube.  Loading  and  unloading  of 
substrate  wafers  is  accomplished  through  a  load-lock  chamber  beneath  the  radio 
frequency  (rf)  section  of  the  reactor  that  can  be  evacuated  by  a  turbomolecular 
pump.  After  purging  with  ultra-pure  nitrogen,  sample  transfer  can  be  executed  us¬ 
ing  a  magnetic  transfer  rod.  Gas  is  purging  through  the  gap  between  the  susceptor 
and  the  reactor’s  base  to  avoid  flow  of  gas  mixtures  to  the  mechanical  workings 
behind  the  susceptor.  The  quartz  glass  reactor  is  connected  at  the  inlet  to  a  source 
vapor/process  gas  flow  control  and  switching  panel  that  directs  individual  streams 
of  source  vapor  saturated  carrier  gas  either  to  a  vent  line  or  to  the  reactor.  Thus, 
pulsed  operation  separating  plugs  of  source  vapor  saturated  carrier  gas  by  plugs  of 
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high  purity  carrier  gas,  flow  rate  modulated  flow  or  continuous  flow  can  be  imple¬ 
mented  for  all  source  vapors  without  change  in  reactor  pressure  or  total  flow.  Two 
optical  windows  at  the  Brewster  angle  of  the  substrate  are  attached  to  the  sides  of 
the  reactor.  They  allow  for  the  real-time  process  monitoring  utilizing  p-polarized 
reflectance  spectroscopy  (PRS)  (see  e.g.,  [17]  and  the  references  therein). 


Figure  1 .  Schematic  representation  of  a  horizontal,  quartz  reac¬ 
tor  in  a  steel  confinement  shell 

To  demonstrate  the  feasibility  of  using  POD  technique  as  a  reduced  basis 
method  for  computation  of  feedback  controls  and  compensators  in  a  high  pressure 
CVD  reactor,  we  will  restrict  our  study  to  a  two-dimensional  rectangular  domain 
(Fig.  2)  representing  the  longitudinal  cross  section  through  the  center  of  the  reac¬ 
tor.  We  consider  the  deposition  of  InP  using  pulsed  trimethyl-indium  (TMI)  and 
phosphine  (PH3)  as  source  vapors  and  hydrogen  as  carrier  gas.  In  particular,  at 
first  only  carrier  gas  flows  through  the  reactor.  After  the  flow  reaches  steady  state, 
a  pulse  of  reactant  (e.g.,  TMI)  diluted  with  carrier  gas  enters  the  reactor.  After  the 
pulse,  the  reactor  is  then  flushed  with  carrier  gas.  This  process  is  then  repeated 
for  another  reactant.  Pulsing  of  the  III-V  source  materials  prevents  nucleation  of 
the  film  in  the  gas  phase  and  make  PRS  observation  and  analysis  possible. 

We  consider  only  trace  amounts  of  reactants  mixed  with  the  carrier  gas.  Un¬ 
der  this  dilute  approximation,  we  can  classify  CVD  processes  as  a  quasi-transient 
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Figure  2.  Two-dimensional  cross  section  of  the  CVD  reactor 
(height=0.011m,  length=0.156m,  substrate  length=0.048m) 


flow  (steady-state  flow  with  transient  species).  The  steady-state  flow  is  described 
by  the  following  set  of  equations  [12] 

(continuity) 

V  •  {pv)  =  0,  (1) 

(momentum) 

pv  ■  Vv  =  —VP  +  V  •  r  —  pg ,  (2) 

where  the  viscous  stress  tensor  is  of  the  form 

2  _  _  _ 

r  =  —  -p(V  ■  v)I  +  /i(W+  WT),  (3) 

O 

(energy) 

pcpv-  VT  =  V-(fcVT),  (4) 

where  g  is  the  gravitational  acceleration,  u,  T,  and  P  are  the  velocity,  temperature, 
and  pressure,  p,  cp,  and  k  are  the  viscosity,  specific  heat,  and  thermal  conductivity 
of  the  carrier  gas.  The  density  variations  are  modeled  as  [12] 

p  =  po[l-0(T-To)],  (5) 

where  To  is  a  reference  temperature,  po  is  a  reference  density  calculated  from  the 
ideal  gas  law  at  the  reference  temperature  and  reactor  pressure,  and  f3  is  the  volume 
coefficient  of  expansion  (/3  =  1/T).  In  addition,  we  consider  a  hydrogen  carrier 
gas  at  atmospheric  pressure.  Temperature  dependent  values  for  p,  k,  and  cp  are 
linearly  interpolated  from  measurements  taken  from  the  available  literature  [12]. 
A  parabolic  velocity  flow  profile  is  specified  at  the  inlet  (IT),  with  an  average  inlet 
velocity  of  0.1147  m/s.  No  slip  (zero  velocity)  boundary  conditions  are  imposed 
on  those  portions  of  the  model  corresponding  to  the  reactor  walls  (T2,  T4,  T5,  and 
T6).  Room  temperature  boundary  conditions  are  imposed  at  the  inlet  and  along 
the  upper  wall  (T2).  Along  the  bottom  wall,  the  substrate  (T5)  temperature  is 
fixed  at  800°  K  ,  with  a  non-linear  temperature  decrease  from  the  substrate  edge 
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to  the  inlet  (T6)  and,  similarly,  from  the  substrate  edge  to  the  outlet  (T4)  (see 
Fig.  2). 

Steady-state  solutions  for  v,  T,  and  p  obtained  from  equations  (1-4)  are  then 
used  in  the  time-dependent  species  equations  for  the  precursor  mass  fractions  [12], 

dY  i  /  ,  Nr 

+  v  ■  \7Yn  =  -  V  •  (pDn VYnj  +  £Yn<  ,  (6) 


where  Dn  is  the  diffusivity  of  the  species,  Yn  is  the  mass  fraction  of  the  nth  species, 
Nr  is  the  number  of  gas  phase  reactions,  and  rni  is  the  rate  of  production  of  species 
n  in  the  zth  chemical  reaction. 

Under  the  reactor  conditions  considered  here  (H2  carrier  gas,  800°K  substrate 
temperature,  and  1  atm  pressure),  there  are  no  effective  gas  phase  reaction  mech¬ 
anisms  for  phosphine,  and  the  only  significant  gas  phase  reaction  for  TMI  is  the 
decomposition  of  TMI  to  MonoMethyllndium  (MMI)  and  two  methyl  molecules, 
In(CH3)3  — >  InCEU  +  2CH3.  This  reaction  can  be  described  as  a  first-order  Ar¬ 
rhenius  reaction  [12] 

rn  =  "n^71  k0e(~E/RT)YTMh  (V 

vvtmi 

where  vn  refers  to  the  stoichiometry  of  species  n  in  the  reaction,  Wn  and  ITtmi 
refer  to  the  molecular  weight  of  species  n  and  TMI,  respectively,  ko  =  5.25  x  1015 
s-1  is  the  rate  constant,  and  E  =  47.2  kcal/mol  is  the  activation  energy  [18]. 

The  tracking  control  problem  that  we  formulate  is  to  find  the  mass  fractions 
of  TMI  and  phosphine  at  the  inlet  (Ti)  of  the  reactor  in  order  to  obtain  a  desired 
flux  qr(t)  of  reactants  at  point  xp  on  the  susceptor  (T5).  That  is,  we  consider  to 
minimize  a  cost  functional  of  the  form 


subject  to 


J(u) 


[u'Ru  +  (q  —  qrYQiq 


qT)]  dt , 


+*‘VYn 


Yn(  0,£) 
Yi{t,x) 
Y2(t,x) 
Y3(t,x) 
Y„(t,x) 

dYn  ( t,x ) 
dn 


n 


jV  ■  ( pDnVYn )  +  A nY] 

Vn0(x) 

Ui(t) 

U2{t) 

0 

0 

0 

1,2,3, 


on  n 
on  n 
on  n 
on  T5 

on  T2  U  T3  U  T4  U  T6 


(8) 


(9) 


where  Yi,Y2,  and  Y3  refer  to  the  mass  fractions  of  TMI,  phosphine,  and  MMI, 
respectively;  ui(t),u2(t)  are  the  controls  corresponding  to  TMI  and  phosphine; 
A r(T)  =  -k0e(~E/RT\  A 2(T)  =  0,  and  A 3(T)  =  (W3/W1)k0e^E/RT^,  v,  p,  T  are 
the  steady  state  solutions  to  equations  (1-4);  and  denotes  the  outward  normal 
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derivative.  Finally,  the  general  flux  vector,  q(t),  at  the  point  xp  is  given  by 


q(t) 


Qln(t) 

qp(t ) 


D 


Win 

1  W 1 


§Xl 

dn 


d2 


-  + 

Xp 

Wp  dY2  I 
W2  dn  I 


Wi„ 

W3 

xp 


0Y3 
dn  xp 


(10) 


where  Wjn  and  Wp  are  the  molecular  weights  of  indium  (a  component  of  Y\  and 
I3)  and  phosphorus  (Y2),  respectively. 

We  note  that  since  the  methyl  (CH3)  molecules  do  not  participate  in  film 
growth  or  otherwise  affect  the  transport  properties  (under  the  dilute  approxima¬ 
tion),  we  do  not  include  them  in  the  state  equations  (9).  The  reactor  walls  (T2, 
T4,  and  T6)  are  assumed  non-absorbing,  and  the  substrate  (T5)  is  assumed  to  be 
perfectly  absorbing  (concentration  of  zero) .  Temperature  dependent  values  for  the 
diffusivities  Dn  in  hydrogen  are  linearly  interpolated  from  values  taken  from  the 
available  literature  [18]. 

We  next  use  a  penalty  boundary  formulation  on  the  species  state  equa¬ 
tions  (9)  to  change  all  Dirichlet  boundary  conditions  to  Neumann  conditions.  For 
example,  the  Dirichlet  condition  Y-\  (t,  x)  =  u±(t)  is  reformulated  as  ^  = 

l(Yi  (t,  x)  —  where  e  is  a  small  parameter  (for  most  of  our  calculations  we 

used  e  =  ICC3).  This  boundary  conditions  reformulation  provides  a  natural  setting 
for  the  Galerkin  procedure  as  well  as  for  the  control  formulation.  More  specifically, 
writing  the  state  equation  in  weak  form  using  test  function  Wj,  integrate  by  parts, 
and  applying  the  modified  Newmann  conditions,  we  obtain 

fn~dtLwj  ^  =  —  In  (y  '  wj  dQ  ^  Jq  Dn^7Yn  ■  S7vjj  dfl 

+  Jq  p  wjDnVYn  ■  Vpdfl  +  Jq  A nYnwj  dfl  (H) 

T  f  /|-  |  j  -  djj  rdrl  Y„  ds  —  Jj-  j  /djjDnun  ds  , 

where  n  =  1,2,3  and  U3  =  0.  The  mass  fraction  of  the  nth  species  is  approximated 
as  a  linear  combination  of  the  Mn  most  significant  POD  basis  elements  as 


Mn 

Y™n(t,x)  =  (£)$„,*(£), 

i—1 


(12) 


where  Mn  <C  N  and  is  the  zth  POD  basis  element  corresponding  to  the 
nth  species.  Using  the  representation  (12)  and  the  orthonormality  of  the  {<f>}'s, 
we  apply  a  Galerkin  procedure  to  the  weak  form  (11)  to  obtain  a  system  of  M 
ordinary  differential  equations  for  the  coefficients  yn^ 


yM(t)  =  AMyM(t)  +  BMu(t),  u(t) 


Ul(t) 

u2  (t)  J  ’ 


(13) 


where  M  =  X)n=i  AM  is  an  M  x  M  matrix,  and  BM  is  an  M  x  2  matrix. 

We  remark  that  the  POD  modes  for  each  species  are  constructed  from  150 
snapshots  ( N  =  150)  taken  during  the  three  second  cycle  (2  s  pulsing,  1  s  clearing) 
of  each  source  species.  Each  solution  vector  represents  the  species  mass  fraction 
at  the  453  nodal  points  and  corresponds  to  a  time  increment  of  0.03-s  in  the  time 
range  from  0  to  3  seconds.  However,  for  the  reduced  order  model,  we  used  only 
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19  POD  modes  (M  =  19:  MpMI  =  8,  ^Phospine  =  ^MMI  =  3),  which  yield 
a  worst  captured  variability  of  99.995%  for  MMI  (X^=i  ^MMIf/^fc=*i  ^MMIfc  = 
99.995).  While  the  captured  variability  suggests  the  proper  order  for  accurate  re¬ 
duced  order  model  simulations,  there  are  additional  order  questions  related  to  the 
control  system  to  be  used  in  determining  reduced  order  gains  and  compensators. 
For  example,  the  ranks  of  the  controllability  and  observability  matrices  have  some¬ 
times  been  found  to  be  useful  criteria  to  help  in  the  choice  of  the  number  of  modes 
to  use  in  control  design  applications  for  the  reduced  basis  representation  [11,  12]. 

The  reduced  order  species  state  model  (13)  is  linear  in  both  the  state  and 
control  variables  and,  therefore,  linear  control  methodologies  can  be  applied  to 
obtain  the  optimal  feedback  tracking  control.  However,  the  application  of  this 
feedback  control  to  the  reduced  order  model  requires  a  linear  state  estimator  since 
only  partial  state  observations  of  the  fluxes  of  In  and  P  at  the  substrate  center  are 
available.  These  partial  state  observations  are  compatible  with  current  PRS  sensing 
technology.  For  detailed  calculations  of  the  POD  basis  elements,  the  reduced  order 
model,  and  the  compensator-based  optimal  feedback  tracking  control  see  the  recent 
article  [12]. 

To  show  how  well  the  control/compensator  system  designed  above  via  the 
reduced  order  model  performs  when  used  in  the  actual  physical  experiments,  we 
computationally  test  the  reduced  order  control/compensator  design  on  the  full 
system,  which  is  approximated  by  453  x  3  quadratic  finite  elements.  We  emphasize 
that  the  reduced  order  model  is  formulated  using  only  19  POD  basis  elements, 
a  substantial  order  reduction  from  453  x  3.  Fig.  3  depicts  plots  of  the  observed 
fluxes  as  functions  of  time.  It  clearly  shows  that  the  system  is  able  to  closely  track 
the  time  dependence  of  the  desired  flux  profile  (shown  in  dotted  line)  without 
significant  delays.  It  also  confirms  that  the  ability  of  the  system  to  match  the 
target  flux  is  sensitive  to  the  design  parameter  Q. 


4.  Nonlinear  Tracking  Control  and  State  Estimator 

4.1.  Tracking  Control  for  Nonlinear  Systems 

In  general,  the  mathematical  model  developed  in  §3  has  to  be  linked  with  a  surface 
kinetics  model  describing  the  decomposition  kinetics  of  the  organometallic  precur¬ 
sors  involved  and  their  incorporation  into  the  film  deposition.  For  the  growth 
of  epitaxial  GaP  heterostructures  on  Si(001)  substrates,  a  reduced  order  surface 
kinetics  (ROSK)  model  is  proposed  in  [19].  This  surface  kinetics  model  is  a  non¬ 
linear  system  of  ordinary  differential  equations  and  when  combining  with  the  gas 
phase  model  will  yield  a  system  of  nonlinear  differential  equations.  In  this  section, 
we  summarize  our  development  of  nonlinear  estimators  and  nonlinear  feedback 
tracking  controls  that  are  applicable  to  a  wide  class  of  systems  including  the  high 
pressure  CVD  systems  as  considered  here  [20]. 
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Figure  3.  Observed  fluxes  as  a  function  of  time  for  different 
values  of  the  control  parameter  Q:  small  (solid  line),  medium 
(dashed  line),  and  large  (dash-dot  line).  The  target  flux  profile 
is  also  shown  (dotted  line)  for  reference.  The  system  is  able  to 
closely  track  the  desired  flux  profile. 


Consider  the  following  nonlinear  control  system 

(  x(t)  =  f(x(t))  +  Bu(x(t),t) 

<  x(0)  =  x0 

l  y(t)  =  Hx(t), 

where,  for  this  presentation,  the  tracking  variable  y  is  taking  to  be  a  linear  function 
of  the  state  variables  (see  [20]  for  the  nonlinear  tracking  variable  case).  In  addition 
to  the  nonlinear  state  equation,  the  cost  function  for  the  tracking  problem,  with 
a  desired  trajectory  r(t),  is  given  by 

J(xo,  u)  =  -  J  ((y  -  r)T Q{y  -  r)  +  uT Ru)  dt. 

Now  rewriting  the  nonlinear  function  as  f{x)  =  A(x)x  and  solving  the  necessary 
optimality  conditions,  we  obtain  the  optimal  feedback  control 

u(x,  t )  =  ~R^1Bt  [II  (x(t))  x(t)  +  s(t,  x)] , 


(14) 
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where  n(cr)  is  the  solution  to  the  so-called  state  dependent  Riccati  equation 
(SDRE) 

n  ( x )  A(x)  +  at(x) n  (x)  -  n  (x)  br~1bt n  (x)  +  htqh  =  o  (15) 


and  s(t,x)  is  the  solution  to  the  following  two-point  boundary  value  problem 

8  =  -AT{xnom)s  +  n(®„  'nn-lRT-^ 


)BR~1BTs  +  HTQr 

T 


-  Elil  ( SrnmOi  (fc|’,(^m))  (n  (a 

D-fTl  (^nom)  Xnorn 

A(x  norn^)  X  norn  BR  B  (II  ( xnorn )  xnor 


n)  %ri 

s) 


+  s) 


(16) 


with  xnorn(0)  =  Xq  and  s(Tf ,  xnom(Tf))  =  0,  (see  [20]  for  details). 

The  SDRE  (15)  is  solved  using  a  power  series  approximation.  We  begin  by 
splitting  A  into  a  constant  part  and  a  state-dependent  part  as  A(x)  =  Aq  + 
eg(x)AAc,  where  e  is  a  temporary  variable  used  for  the  expansion  that  will  be  set 
to  1  later.  We  next  write  II  as  a  power  series  in  e,  as 

OO 

n  {x,e)  =  Y,en9n{x){Ln)c,  (17) 

n—0 

where  II  as  well  as  each  ( Ln)c  is  symmetric.  Substituting  these  expansions  into  the 
state-dependent  Riccati  equation  (15)  and  matching  terms  with  the  same  powers 
of  e  we  obtain  the  following  set  of  equations  for  determining  the  const  ant- valued 
matrices  ( Ln)c : 

(■ L0)cA0  +  A^(L0)c-(L0)cBR-1BT(L0)c  +  Q  =  0(18) 
(ir)c  (A)  -  BR~1BT(L0)c)  +  (Aq  -  (L0)CBR~1BT)  {L{]c 

+(L0)cAAc  +  AAq(L0)c  =  0  (19) 
{Ln)c  (^o  -  BR~1BT(L0)c )  +  (Aq  -  (L0)CBR~1BT)  (Ln)c 

n—  1 

+(Ln-i)cAAc  +  AATc(Ln_1)c-YJ(,(Lk)cBR-1BT(Ln_k)c )  =  0.(20) 

fe= l 

Equation  (18)  is  the  standard  Riccati  equation  for  the  linear  part  of  the  system,  Ao, 
which  can  be  solved  easily.  Equations  (19)  and  (20)  are  constant-valued  matrix 
Lyapunov  equations,  for  which  stable  and  efficient  algorithms  also  exist  in  the 
literature. 


4.2.  State  Estimation  for  Nonlinear  Systems 

We  consider  a  nonlinear  control  system  with  a  nonlinear  measurement  of  the  form 
f  x(t)  =  f(x(t))  +  Bu(xe(t),t) 

\  z(t)  =  c(x(t)). 

The  control  for  a  tracking  problem  is  given  by 

u(xe,  t)  —  R  B  (II  (xe)  Xe  “t“  s(tj  Xriom')') 


as  discussed  in  the  last  section  except  now  in  terms  of  the  estimated  state  xe. 
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The  estimated  state  will  be  formulated  by  an  ordinary  differential  equa¬ 
tion  similar  to  the  state  equation,  with  a  gain  matrix  (found  using  a  dual  state- 
dependent  Riccati  equation)  applied  to  the  difference  between  the  measurements 
of  the  actual  and  estimated  states.  The  coupled  actual  and  estimated  states  are 
given  by 

(  x  =  A{x)x  -  BR~lBT  [II(a;e)a;e  +  s(t,  xnom )] 

<  xe  =  A(xe)xe-  BR~1BT  [U(xe)xe  +  s(t,Xnom)}  (21) 

[  +  L(xe )  [2  -  c(xe)\ , 

where  the  state  estimation  gain  is  found  by 

L(xe)  =  E(a ^(CofV-1  (22) 

with  E(xe)  is  the  solution  to  the  dual  state-dependent  Riccati  equation 

Z(xe)AT(xe)  +  A{xe) E(xe)  -  E(cce)(C,o)TV'_1C'oE(a;e)  +  U  =  0.  (23) 

For  the  purposes  of  finding  the  estimator  gain  in  equations  (22)-(23)  the  nonlinear 
measurement  function  z(t)  =  c(x(t))  is  rewritten  as  matrix  function  multiplication 
c(x)  =  C(x)x  and  to  choose  Co  =  C(0).  We  note  that  the  nonlinearity  of  the 
measurement  function  does  remain  in  the  estimator  system  (21)  itself,  and  the 
nonlinearity  of  the  system  dynamics  remains  in  (21)-(23).  The  estimator  gain 
SDRE  (23)  is  solved  using  the  power  series  approximation  in  an  analogous  fashion 
as  described  in  previous  section. 

4.3.  A  Flight  Dynamics  Example 

We  apply  the  nonlinear  estimator  and  feedback  tracking  control  presented  in  pre¬ 
vious  sections  to  the  flight  dynamics  example  from  [21].  The  control  system  is 
given  by 

x  =  (A0  +  x2Anl)  x  +  Bu,  z(t)  =  [#i,  X2,x5]t 


where  the  matrices  Aq,  A^l 

and  B  are 

given  by: 

'  -0.0443 

1.1280 

0.0 

-0.0981  0.0 

-0.0490 

-2.5390 

1.0 

0.0 

-0.0854 

A0  — 

-0.0730 

19.3200 

-2.2700 

0.0 

22.6834 

0.0490 

2.5390 

0.0 

0.0 

0.0854 

0.0 

0.0 

0.0 

0.0 

20.0 

'  -0.2317 

0.0 

0.0 

0.0 

0.0 

-1.2760 

-0.7922 

0.0 

0.0 

0.0206 

Anl  = 

0.1020 

64.2940 

-13.9710 

0.0 

-5.4167 

1.2760 

0.7922 

0.0 

0.0 

-0.0206 

0.0 

0.0 

0.0 

0.0 

0.0 

B  =  [ 

0.0  0.0 

0.0  0.0 

20.0  ]T . 

The  cost  functional  to  be  minimized  is 

B 

o 

II 

tO  1  H 

/  ((»- 
Jo 

r)T Q(y  -  r 

)  +  utRu )  dt, 
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The  state  variables  in  this  model  represent  the  flight  conditions  of  the  aircraft:  x\ 
is  the  deviation  of  the  velocity  from  the  level  flight  trim  value  of  l(100m/s)  (given 
in  units  of  (100m/s)),  X2  is  the  deviation  of  the  angle  of  attack  from  the  trim  value 
of  4.2(7t/180)  radians,  X3  is  the  pitch  rate  in  rad/s,  £4  is  the  flight  path  angle  in 
radians,  and  x§  is  the  deviation  of  the  canard  deflection  angle  in  radians  from 
the  trim  value,  which  is  not  given.  The  control  u  is  the  input  canard  deflection 
in  radians.  The  canards  are  control  flaps  which  can  deflect  downward  by  up  to 
90(7t/180)  radians.  Finally,  y(t)  =  X4 (t)  is  the  tracking  flight  path  angle  and  r(t) 
is  the  desired  flight  path  angle. 

The  weights  are  Q  =  1,  R  =  1,  U  =  IOO/5  and  V  =  I3.  The  actual  state 
starts  at  the  origin,  but  the  estimated  state  starts  slightly  off  the  actual,  at  (xe)o  = 
(0,  0,  0,  5(7t/180),  0)t.  Figure  4  depicts  the  estimated  state  almost  converging  to 
the  actual  state  by  the  time  of  the  desired  X4  increase,  and  remaining  close  to 
the  actual  state  for  the  rest  of  the  time  period.  In  Figure  5  we  plot  the  actual 


Figure  4.  Actual  and  estimated  states  for  nonlinear  tracking 
control/state  estimator 


state  when  controlled  using  our  fully  nonlinear  algorithm,  as  well  as  when  using 
the  linear  SE  gain  control  (as  proposed  in  [22]),  and  the  fully  linearized  control.  It 
can  be  seen  that  the  linear  control  overshoots  significantly  at  the  top  of  the  ascent 
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Figure  5.  Comparison  of  tracking  controls/state  estimators  with 
inaccurate  (xe)o- 

and  is  very  slow  to  return  to  0.  The  other  two  methods  produce  virtually  identical 
results  (the  difference  is  indiscernable  in  the  plots  in  Figure  5). 

In  the  recent  manuscript  [23],  we  have  successfully  applied  these  methodolo¬ 
gies  to  feedback  tracking  control  of  the  GaP  film  thickness  in  a  high-pressure  CVD 
reactor  that  is  more  complex  and  realistic  than  the  one  presented  here.  Some  of 
the  complexities  that  were  considered  in  [23]  include  3-dimensional  flow  region,  10 
atm  operating  condition,  multiple  species  that  include  gas-phase  kinetics  as  well 
as  nonlinear  surface  kinetics,  and  nonlinear  partial  state  observation. 

5.  Conclusions 

In  this  paper  we  report  on  the  development  of  nonlinear  compensators  and  non¬ 
linear  feedback  tracking  control  methodologies  that  can  be  applied  to  high  pres¬ 
sure  CVD  systems.  We  also  present  successful  computational  implementation  of 
reduced  order  feedback  control  of  pulsed  high  pressure  CVD  III-V  film  growth 
involving  the  transport  of  multiple  species  with  linear  gas  phase  reactions.  The 
combination  of  reduced  order  model  methods  based  on  proper  orthogonal  decom¬ 
position  techniques  with  nonlinear  compensator-based  feedback  tracking  control 
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can  provide  a  powerful  tool  for  treating  more  complex  situations  that  may  be  en¬ 
countered  when  nonlinear  gas  phase  and/or  nonlinear  surface  phase  reactions  are 
present. 
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